A physics-based model of swarming jellyfish

We propose a model for the structure formation of jellyfish swimming based on active Brownian particles. We address the phenomena of counter-current swimming, avoidance of turbulent flow regions and foraging. We motivate corresponding mechanisms from observations of jellyfish swarming reported in the literature and incorporate them into the generic modelling framework. The model characteristics is tested in three paradigmatic flow environments.


Introduction
Scyphozoans populate the oceans since the late proterozoic eon and according to geological records they are one of the first multicellular organisms on planet Earth [1]. The life cycle of scyphozoans is complex including in many species the tiny, cryptic, benthic stage of the polyp. But, most prominent is the planktonic stage of the medusa, better known as jellyfish [2,3]. Despite their appearance, medusae are among the most efficient swimmers in the oceans [4,5] and their massive occurrence in different habitats has major impacts on the biosphere, tourism, and economics [6][7][8][9]. While jellyfish have existed throughout human civilization, the adverse effects of climate change on the ocean habitats are hypothesized to intensify jellyfish blooms and the associated impacts [6,10].
The fact that scyphomedusa establish swarms and interact with humans on numerous levels has motivated researchers to investigate the relationships between large jellyfish aggregations and environmental conditions [11][12][13][14][15][16]. Results of these studies have mainly contributed to empiric models for prediction of jellyfish proliferation in estuaries and bays [17][18][19] or along open coasts [20][21][22]. These approaches represent the coarse end of a spatio-temporal continuum of models, ranging from regional scale to the level of single-agent swimming. Indeed, many studies have examined the efficiency of different propulsion mechanisms in single medusae, experimentally [23][24][25] and computationally [26][27][28][29][30][31][32][33], based on detailed observations and direct numerical simulations of bell oscillation and the resulting vortex dynamics in the surrounding water.
In addition, jellyfish respond to a variety of environmental stimuli, affecting the emergence and movement of swarms as a whole. Observations suggest that jellyfish modify their swimming behavior with respect to: the location of other individuals [34], the water temperature [13,18], salinity [11,14], advection by water currents, turbulence [19,23]  prey [35,36]. However, there is a large mismatch, between what is known regarding jellyfish ecology and the accuracy of jellyfish swarm prediction based on circulation and distribution models. Thus, there is a great need for a theoretical modelling framework that ties together single-agent dynamics, response mechanisms and swarm behavior, applicable in large scale circulation models [3,[37][38][39][40]. Thus, in this paper, we pursue two research goals: First, we propose a mechanistic model framework for the swarming dynamics of jellyfish based on active Brownian particles (ABPs) [41,42]. A key aspect of the model is that it reaches far beyond existing approaches for modelling of jellyfish proliferation (see Table 1 below): Jellyfish in the model are considered to move actively and they can make decisions based on environmental information. In particular, we pursue the idea that the environmental stimuli are processed by the neuronal network of a jellyfish. This network generates a coherent strong pulsation rhythm which gives rise to the bell movement in many jellyfish species. Accordingly, we consider a jellyfish in our model to be a moving nonlinear oscillator, or swarmalators [43][44][45]. A swarmalator is the next step towards modelling of agent-based swarm phenomena, combining active-matter research and the theory of nonlinear oscillators [46][47][48][49]. In particular the latter allows us to incorporate several important biological features into the model, for example the species-dependent response to external stimuli.
The second and main research goal is to adapt the introduced parameters of the model according to the existing literature and experience from field observations [23,35,65]. On the one hand, the particle-based description allows to extract a consistent mathematical description of individual behavior by means of data analysis techniques for oscillatory systems [66][67][68][69][70][71]. On the other hand, the ensemble dynamics can be tested in a transient development process of prediction and continuous parameter adjustment, using physics-informed machine learning [51,[72][73][74] where the parameter set of this work serves as a starting point. Here we integrated data collected for species of the orders Rhizostomeae and Semaeostomeae, species which are known for their swarming behavior. Despite the differences in the bell shape (prolate vs oblate) and swimming patterns (jet-propulsion vs rowing-propulsion) in adult medusae of the two orders, we have managed to provide an overarching set of parameters. However, given the fact that the fields of active matter dynamics and jellyfish ecology have been largely disconnected up to now, a certain amount of intuition is needed to construct an overarching model and to design numerical tests for its validation. As a consequence, we exploit only a minimal set of mechanisms, related to environmental inputs, while we leave out most of the agent-agent interactions. Nevertheless, the resulting active oscillatory swimmers are able to mimic biological phenomena found in swimming jellyfish [75]. Our analysis and characterization of the model performance focuses on theoretical and computational strategies to investigate pattern formation in networks of moving oscillators from the perspective of active-matter research [52,55,56,[76][77][78][79][80][81] and tries to relate these outcomes to the biological background. We elaborate in the concluion section which theoretic quantities should be extracted from experimental data to strengthen the observational evidence for active-jellyfish modelling.

Methods
We obtain the flow field U(x, t) (hereafter x denotes position on the horizontal plane and t denotes time) from the incompressible Navier-Stokes equations which we simulate by means of a second-order in time and space method on a staggered grid using a Successive Over-Relaxation (SOR) pressure-solver [82,83]. The dynamics of the prey concentration F(x, t) is modelled by an advection-diffusion equation coupled to the flow. We fix the diffusion coefficient of the prey to be D F = 0.001 m 2 s −1 [84].
The trajectories of the jellyfish are simulated by the second-order stochastic Heun-method [85]. We couple the dynamics of the jellyfish to the flow by using a bi-linear interpolation of the local field quantities onto the position of each single jellyfish agent j [86] so that for instance F j = F(x j , t).
For the statistical analysis of the swarming behavior we make use of ensemble averages, the Pearson correlation [87] and the indication number: Here j denotes the serial number of a jellyfish agent and M � N is the overall number of agents considered for averaging. W can be any simulated field variable, e.g. the concentration of prey. N is the total number of agents. In reality jellyfish swarms are composed of hundreds of thousands of individuals in the open sea, or just several few individuals in controlled tank experiments [12,25,88,89]. In this work we use N = 128 agents as a tradeoff between theoretical demands and experimental limitations.
To obtain statistically significant results, we average values of Eq (1) over K = 16 model runs using a Gaussian kernel with a window of 15 seconds [90]. The resulting double average is denoted as hWi � hW M i K . σ 1 and σ 2 are the standard deviations of variables W 1 and W 2 . N ðWÞ measures the fraction of jellyfish that can be found in domains fulfilling the condition ensemble average for passive tracersŴ � W passive;N if not stated otherwise. For example regarding prey searching, we expect that jellyfish will try to maximize their access to ambient prey while passive tracers will not show any response. Thus, when for a majority of jellyfish F j >F holds, N ðFÞ � 1 indicates an excess of preying performance and most jellyfish are situated in regions where the prey concentration is at least higher thanF. An exemplary separation of the flow domain is depicted in Fig 6. In technical terms, the separation is represented by the Heaviside function H which returns unity if W >Ŵ and zero otherwise.

Towards a description of active swarming jellyfish
We model a single jellyfish agent as an active over-damped particle at horizontal position x = (x, y) > and time t, having a velocity v(x, t) and orientation θ(x, t) [41,42]. In the following, we propose a dynamic model for the position and the orientation of the agent which allows us to capture paradigmatic behaviors of jellyfish.
Jellyfish process environmental information by a neuronal network of several thousand neurons. The dynamics of this network can be regarded as a perturbed internal dynamics [91][92][93][94] in which neurons fire coherently and act as a single large oscillator to drive the bell pulsation. Moreover, it has been shown that the neuronal network of jellyfish features properties of a circadian oscillator [89,95]. Thus, the bell pulsation can be regarded as the resulting average network oscillation [5,23,46]. For such a process it has been shown theoretically [46,[96][97][98][99][100][101] and empirically [102,103] that a low-dimensional description, in terms of a phase variable φ(t) is possible and beneficial for a swift but reliable model development. Moreover, phase dynamics models are computationally light, but at the same time a reliable simplification, as they are based on observational evidence and a fit to the theory [66,69,99,103].
Another aspect of jellyfish motion is that swimming patterns can change abruptly in response to stimuli [65,75]. In such cases jellyfish will switch from a state of relative inactivity into a state where the frequency and the strength of the bell oscillation increase and the swimming trajectories encompass a significantly larger volume of fluid [36, 104,105]. Such changes in swimming patterns can even cause a swarm to cross flow barriers [35]. Following such periods of elevated activity, jellyfish will return to a state of relative quiescence. To capture this behavior, we introduce an additional activity variable A which causes parametric switching [101,106,107]. In accordance with theoretical studies on relaxation oscillators, this variable can be regarded as the leading-order amplitude perturbation from a limit cycle [98,108,109]. We assume that this degree of freedom in the jellyfish dynamics is stimulated by external inputs and that it decays exponentially at a rate l A , in order to ensure it returns to a quiescent state.
Following these generic considerations, the dynamics of a single jellyfish j in a swarm of N agents is described by the four differential equations: Here, we have used the ensemble notation X = [x 1 , . . ., x N ], φ = [φ 1 , . . ., φ N ], θ = [θ 1 , . . ., θ N ]. F indicates the prey concentration and U denotes the water horizontal current vector. C ¼ r � U is the vorticity of the water flow. We take its absolute value as a simplified measure for the amount of turbulence in the flow [110]. In the following we exemplify the generic set of Eq (2). For convenience, a scheme of this jellyfish agent model is illustrated in Fig 1. A simplistic comparison of the model to already existing research on active matter and jellyfish dynamics is presented in Table 1. Next we explain the bottom-up rational of the model. Agent activity A. Here we take into account only the presence of prey as an incentive to activate the jellyfish dynamics and assume that the activity evolves according to: The decay parameter l A corresponds to the time for which an agent is active. We fix this value to l A ¼ 0:005 s −1 , equivalent to a characteristic activity time of roughly 3 minutes [65,111]. We assume that the excitation is a linear function of the magnitude of the local prey concentration gradient |rF j |. Additionally, the strength of the linear forcing is modulated by the relative orientation of an agent with respect to the direction of the prey concentration gradient. Defining the angle of that gradient as δ F (satisfying tan(δ F ) = @ y F/@ x F), the term in the square bracket of Eq (3) ensures that only swimming towards the prey (|θ j − δ F | < π/2) evokes an increase in activity. We choose this approach to mimic a basic type of energy saving strategy, by favoring directions towards prey over others.
Bell oscillation φ. The oscillatory movement of the bell is modelled by a phase oscillator [44, 45,48,60,112] with a frequency o j ðA j Þ. Additionally, the oscillator reacts to perturbations comming from other jellyfish according to a phase coupling function H j (.). Included model mechanisms are advection, swimming against the flow, soft-core repulsion, avoidance of walls, bell oscillations, avoidance of turbulence, local random motion and behavioral shifts due to activity changes (presence of prey). External drivers of the model are the fluid velocity U j , the absolute vorticity jCj j the gradient of prey rF j and the prey concentration F j . https://doi.org/10.1371/journal.pone.0288378.g001

PLOS ONE
We allow the natural frequency of the bell oscillations to vary according to the activity level A. In order to mimic this coupling, we use the response function: Here, a and b are the generic parameters of the response to be set later, and S is a stimulus that evokes the response. Generally, Rð:Þ resembles a nonlinear activation function, commonly associated with dynamics in large neuronal networks and employed in machine learning [113]. The reason for this widespread use is that Rð:Þ mimics a quasi-linear response at small stimuli (R � a b S) and a saturation behavior at large stimuli (R � a), typical for physiological limitations in biological systems (see Fig 3 panel (a)). The resulting natural frequency response is defined as: where ω j,0 is the natural frequency of an agent in an inactive state. In this study we fix ω j,0 = ω 0 = 1.2 rad s −1 [23,111] for all agents. Since the bell oscillation frequency in several species can depend on the size of jellyfish, future experimental work is needed to determine the evolution of its distribution for different species during different seasons [50]. Moreover, we seek to allow for a maximal physiological frequency response which we expect to observe when jellyfish escapes its predators. Thus, we fix the frequency response parameter to f φ = 0.75 [114]. We use a single type of response function, Eq (5), throughout this study both, due to the lack of experimental data and in order to simplify our model. We expect however, that based on future observations, a variety of response functions will be required, similar to paradigmatic models of computational neuroscience [106,107]. The bell oscillation is subject to multiple external stimuli, foremost the turbulence of the surrounding water due to the bell strokes of other agents. The resulting perturbations exert stresses on the bell tissue and generally impact the oscillatory states of the bell muscles and with these, the neuronal network, causing a potential shortening or lengthening of the bell oscillation period. This variability in periodicity is captured by the phase coupling function H j (.) in Eq (2) [46,115]. Due to the lack of observations and for the sake of simplicity we set H j (.)�0, assuming no direct response of the bell to the presence of other agents. Future experiments have to be designed to check for such phase couplings and to extract any potential phase coupling from data [66]. In Eq (2), we indicate the intricate nature of coupling by dependencies on X and θ as the strength of coupling will depend on the position of individuals, the pairwise distance and their orientations.
Orientation dynamics θ. Jellyfish tend to orient themselves according to different stimuli but they are also subject to the eddy motion induced by their own swimming. The resulting turbulent eddies cause a certain angular drift L which acts to reorient the jellyfish, so that G j (.) is the angular coupling function. The angular diffusion according to L j is a frequently observed phenomenon on the microscopic [116,117] and even mesoscopic [118] scale.
Although the properties of this diffusion have not been investigated yet for jellyfish, we expect a stochastic approximation of local eddies to be valid. This is because the reorientation is realized by eddies that remain, at least temporarily, attached to a jellyfish [25, 26] before they separate and may lead to communication among different individuals.
Here, we model the phenomenon by a colored noise [119] which allows us to incorporate a simplistic version of correlations and inertial effects of rotation: The parameters λ θ is a measure for the correlation time of the noise and D is the strength of the white noise η j (t). Both, λ θ and D determine the strength of the angular diffusion. We assume that the resulting motion becomes more prominent when jellyfish search for prey, as individuals start to shuffle water towards their oral arms for feeding purposes [105,114,120]. Thus, we allow for an increase of the noise intensity similar to the natural frequencies of agents: For simplicity we assume equal response parameters f θ = f φ and we fix the essentially free parameter D 0 = 0.1rad 2 s −3 , which determines the angular diffusion in an inactive state The angular coupling function G(.) acts like a potential in the noisy dynamics and causes certain directions to be more favorable than others. In this study, we consider three paradigmatic orientational inputs [75]: First, the agents tend to orient themselves such that they counteract the flow in which they are immersed. Whereas the physiological causes of this behavior are still hypothesized, it is a matter of fact that swimming against the currents increases the survival rate [19]. For example, currents may transport jellyfish to undesired regions such as shore lines and swimming against the flow also increases the chances of catching prey. Second, jellyfish try to avoid regions of turbulent flow in order to optimize their swimming performances and to escape shear flows that may harm their fragile body. Third, jellyfish orient towards the prey concentration, in particular when under starvation [35]. In order to model the interplay of these behaviors, we introduce three separate angular coupling functions g j ðy j À d ðU;C;FÞ Þ with respect to the three local directions d ðU;C;FÞ , of current, absolute vorticity gradient and prey: We assume that the respective angular coupling functions possess a single dominant maximum that can be well approximated by a sinusoidal harmonics [97,112,115,121], resulting in the overall angular dynamics: Constructed in this way (for ε ðU;C;FÞ > 0Þ, the interplay of external orientation inputs tries to establish a stable orientation in parallel of the tumbling.
We can check this for example for the simplified angular dynamics in which we assume that tumbling is absent (L j � 0) and that jellyfish orient just due to the flow of the water. In this case, there exist exactly two fixed points y ? Fig 2). The first point corresponds to swimming in the direction of flow and is by construction unstable. This means that small external perturbation of orientation will grow in this point such that the jellyfish rotate out of alignment. The second point however, is stable and corresponds to swimming against the local direction of flow. For the remaining two orientational inputs, a similar rational applies yielding y ?
, corresponding to avoidance of turbulence and a directed swimming towards prey respectively.
The coupling function Eq (11) gives rise to exactly one stable orientation y ? j and its value depends on the coupling parameters ε ðU;C;FÞ and orientations d ðU;C;FÞ . For example, when the coupling function G(.) is dominated by coupling to prey (ε ðU;CÞ � ε F ), agents mainly orient towards prey such that y ? Fig 2). However, it is conceivable that highorder coupling terms are present in the angular dynamics as well. These terms give rise to several permissible orientations of swimming. For example, we can replace gðy À d C Þ ¼ sinð2ðy À d C ÞÞ in Eq (11) to obtain a high-order coupling function (denoted by G h (.), see panel (b) Fig 2). When noise is present (L 6 ¼ 0), jellyfish will be perturbed out out their stable orientations. In that case, it depends on the amplitude of the coupling constants how likely it is to stay in the direction of a stable equilibrium point. If the dynamics is multi stable, agents can switch in between different equilibria [122]. We exploit this fact to introduce directional decision making by allowing for parametric switches: Here, ε ðU;C;F;0Þ are the respective angular coupling parameters in inactive (U; jCj) or active (F) state. In an inactive state, agents orient against the flow and avoid turbulent regions. An increase of activity will initiate the following cascade of events: • The pulsation frequency o j ðA j Þ increases alongside the amplitude of bell strokes and the coefficient of angular diffusion DðA j Þ, causing an increase in the speed oscillations of swimming and a more frequent change of the swimming direction.
• Jellyfish lose interest in turbulence avoidance. Instead, they initially favor orientations towards high concentrations of prey.
• In response to increase in the instantaneous prey concentration, individuals gradually lose their orientation against the flow.
• When a jellyfish reaches high levels of activity, it again loses its interest in swimming towards prey and establishes a state of free floating.
Additional influence on the orientation can arise due to the presence of other agents, resulting in further structure formation according to more interaction terms g j ðx k ; x j ; y k ; y j ; A j Þ in Eq (11) that depend on the positions, orientations and distance of respective agents j and k [46,59,123,124]. To the best of our knowledge, there are no quantitative estimates for a directional interaction function in jellyfish, hence for now we do not incorporate these effects in the model. We keep in mind however that in dense swarms, probably mostly during the reproduction periods, the mutual interaction between the agents may play a dominant role.
Positional dynamics. We consider an over-damped positional dynamics of the form: The first two terms on the right hand side account for advection and active self-propelled motion, whereas the summations account for the swarm-internal velocity component due to the attraction and repulsion between jellyfish agents. The last term addresses external forcing mainly due to the effect of confinements or obstacles.
We model the active velocity of a jellyfish by: Here,ê ¼ ½cosðy j Þ; sinðy j Þ� is the orientation vector of agent j, V 0 is the maximal propulsion speed of the jellyfish in their inactive state and parameters V a , V b characterize speed adaptions due to changes of current and activity. We parameterize the bell oscillation explicitly. Our choice is inspired by the results found in [4,26,30,[125][126][127] and many other publications on jellyfish swimming which all report or assume pulsed swimming based on varying types of bell movement. In this paper, we use a simplified symmetric version of bell movement that can be adjusted by shape parameter J (see Fig 3(b)). When J is small the self propulsion never fully vanishes, as is the case for species that benefit from inertial effects. In contrast, when J is large, jellyfish temporarily come to a halt as can be observed in species that have developed certain types of jet propulsion. Attraction and repulsion of agents are taken into account by the terms I attr (.) and I rep (.), where the dependencies on θ and φ allow to incorporate physiological abilities of single medusae to sense and communicate in a swarm and in a complex environment. Nonetheless, due to the lack of experimental data on these dependencies we keep both I attr (.) and I rep (.) to be only functions of relative distances. Moreover, since observational evidence for attraction of agents is limited, we assume that the active agents are soft spheres and restrict the interaction to within a radius of R i (see Fig 3(c)) [47,57,77,123,128,129]. Accordingly, N j is the number of agents found within an interaction distance of R i . The resulting repulsion corresponds to unidirectional sensing and communication: We denote the resulting summation of contributions in Eq (14) by I(x j ).
The type of interaction represents a fluid-like repulsion while studies on active colloids often consider molecular interactions (for instance Lenard-Jones potentials) [77,[129][130][131]. We choose this approach because jellyfish are soft and most of the time avoid bumping into each other such that the interaction involves only the fluid. As such, the given interaction resembles the leading order terms of a flow field around an active swimmer on the micro scale [64,132]. Thus, the repulsion I(x j ) only partly represents the turbulent fluid dynamics of jellyfish swimming and we rather see the coupling as a computationally less costly parameterization of the true physical interaction, similar to [43].
Finally, the only external force F ext accounts for interactions with obstacles. In our case, we mimic the avoidance of walls [75,133] in a tank. We model the wall forces by a one-dimensional soft-sphere repulsion, similar to Eq (16) with a repulsion radius of 0.1 m. A simplistic comparison of the overall model and already existing research on active matter and jellyfish prediction models is shown in Table 1.

Discussion and results
We consider three paradigmatic flow environments in a rectangular tank that allow us to compare the performance of our jellyfish model with observational findings. In the first setting, we simulate a cavity flow driven by the movement of the upper wall boundary from the left to the right (see panel (c) of Fig 4). We fix the wall velocity to 0.4 m s −1 [134]. The second setting is a channel flow in a tank of 5 m in width and 40 m in length (see panel (f) of Fig 5). The fluid enters the domain at +20 m and leaves the domain at −20 m. We assume a fully developed Poiseuille flow at the inlet with a magnitude of −0.045 m s −1 in accordance with experimental parameters in [23]. In the third setting we simulate a cavity flow in a tank of 10 m width and 5 m length. We let the left and right wall of the tank move upwards at a speed of 0.4 m s −1 . This way we generate two counter rotating main gyres that separate the domain dynamically into two parts. Throughout all of our simulations we use fluid parameters given in Table 2. Thus, the Reynolds numbers of the cavity flows are Re = 10 5 and for the channel flow it is Re = 4500.

Wall effects
Here, we first discuss what effects arise due to repulsive interactions at walls (F ext ) and due to agent-agent interaction (Iðx j Þ in Eq (14). For this, we simulate two ensembles with reduced

PLOS ONE
The first ensemble mimics passive tracers. They get advected by their local fluid velocities U j and they are confined by the wall forces F ext . In the second ensemble, also volume exclusion is present due to I(x j ). It is well known that in the presence of volume-excluding agent-agent interactions, states of (quasi-) periodic spacial order may emerge [57,135,136]. We seek to avoid such states as they would mask the structure formation process due to swimming and orientation in our simulations. In our model, the transition towards ordered states is influenced by the interaction radius R i and by the number of agents in the domain, which we have to select accordingly.
First, we set the interaction distance at walls for this experiment to 0.5 m. Example results are depicted in Fig 4. We calculate the correlations of velocity and the hexagonal order parameter for the active and passive tracers. To obtain the local hexagonal order, we find for each agent j its n nearest neighbors. Then, δ is the relative angle between agents j and k (obeying tan(δ j.k ) = (y j − y k )/(x j − x k )) and the complex order parameter is given by an average of exp(i6δ j,k ) over these n nearest neighbors. To obtain the ensemble order, we finally average over the ensemble of N agents [136,137]. The velocity correlation measures how much the ensemble follows the background flow while the hexagonal order parameter quantifies local lateral order. For

PLOS ONE
example if the agents were to assemble in perfect hexagons, angles of the 6 nearest neighbors δ j,k = π/6 inside the domain. Accordingly, Hex 6 (X) = 1. If any other symmetry would arise or any other disordered state, it would be indicated by a certain lower level of Hex 6 (X). We see that the swarm undergoes a phase transition from almost free floating to a disordered crystal phase (panels (a) and (b) in Fig 4), depending on the parameters R i and N. This can be observed in Fig 4(d) and 4(e) for the time averages of ΔCorr U,V = hhCorr U,V − Corr U,V;passive ii t and ΔHex 6 (X) = hhHex 6 (X) − Hex passive,6 (X)ii t after an initial transient of 50 s.
At small swarm sizes the agents essentially follow inertial trajectories. When the number of agents increases, we see that for larger interaction radii the velocity correlation drops due to swarm-internal repulsive movements. Depending on R i , the particles essentially lose their ability to float freely. However, at sufficiently large area densities, agents are so closely packed that local repulsive interactions cancel each other (see S1 Movie). Since the positional dynamics Eq (17) is over-damped, it means that particles are still tied to their local fluid velocities such that their motion gains back its correlation with the local flow [135].
In parallel, hexagonal order reveals that for R i = 0.1 m, passive and repulsive swarm behave essentially similarly. At high area densities, we observe that repulsive movement causes a decrease in relative order. On the contrary, we see that the order shows a maximum before it decreases again for larger repulsion radii. This is because the natural increase in hexatic symmetry is disrupted by the flow perturbations and by the quadratic symmetry of the boundary which infiltrates into the domain once long range correlations establish due to dense packing.
Based on this discussion, we set the number of agents to N = 128 and the repulsion radius to R i = 0.1 m as these parameters guarantee almost free floating of the jellyfish-like particles.

Swimming against the flow
We are in particular interested in the swimming dynamics of the jellyfish Rhopilema nomadica. For this species, the ability to counteract the underlying current has been measured in a long channel at different inflow velocities [23]. The swimming speed in the absence of a flow was estimated to 0.067 m s −1 and an active swimming speed of 0.083 m s −1 was reported for an inlet velocity of 0.045 m s −1 .
We adapt the model of passive-repulsive tracers Eq (17) by allowing for an angular dynamics affected by angular diffusion and orientation against the flow in Eq (11), a non-zero swimming speed according to Eq (15) and bell pulsation Eq (6), resulting in the model: Panels (a-e) Fig 5 show snapshots of the resulting swarm dynamics. S2 Movie depicts the time evolution of the swarm.
The, swimming dynamics is defined by the parameters λ θ (inverse of the correlation time for the angular noise), D 0 (amplitude of the angular white noise), V 0 (amplitude of swimming speed in the inactive state), V a , V b (velocity response parameters), ε U,0 (angular coupling constant) and J (shape parameter of velocity oscillations). Here we adjust parameters of the model to experimental results [23]: We choose J = 1 to mimic inertia of medusae during a bell stroke  The corresponding averaged swimming velocity is 0.087 m s −1 . We compensate the excess by allowing for angular diffusion in the simulations. Then, the dynamics depends on parameters ε U,0 , λ θ and D 0 . We fix D 0 = 0.1 rad 2 s −3 and let λ θ 2 [0.2, 1, 5] s −1 (s L ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi D 0 =l y p 2 ½0:71; 0:32; 0:14� rad s −1 ). We perform a parameter scan for ε U,0 2 [0.005�2 (0, 1, . . .7) , 1] rad s −1 and compute the time average of jh _ xi À Uj in the central two third of the simulation time (see panel (g) Fig 5). We divide the average by 0.067 m s −1 to obtain the fraction of swimming speed in presence of a current, similar to the analysis in [23]. Panel (h) in Fig 5 depicts the transition from noisy advection to active counter-current swimming. Our agents swim 21% faster against the flow at ε U,0 = 0.16 rad s −1 and for λ θ = 5 s −1 (Fig 5 bigger triangle in panel h).
The essence of the counter-current swimming can be understood by means of the angular dynamics in Eq (19). We see that the noise has a standard deviation of σ L such that we expect the existence of the fixed point solution y ? j À d u ¼ p to become more likely when ε U,0 > σ L . For the given parameter, this critical average coupling is σ L = 0.14 rad s −1 which indeed is slightly beyond the transition zone of panel (h) Fig 5. In turn, phase slips across the fixed point solution when the noise is temporarily too strong, can still occur. Then, the agents start to tumble or even rotate for some time such that their effective propulsion speed against the flow is drastically reduced (see S3 and S4 Movies). From panel (h) Fig 5 it can be seen that such events are unlikely before the fixed point threshold is reached. Thus, we consider the reconstructed set of parameters to be optimal as it stipulates the frequently encountered phenomenon of criticality in neuronal responses and biological systems in general [138,139].

Turbulence avoidance and structure
The avoidance of turbulent regions in a flow can be implemented in a straight forward way by a second alignment term in Eq (19) of strength ε C;0 : This second term causes the agents to also orient against the gradient of absolute vorticity and thus to swim away from regions of turbulence. As a consequence, flow structures are imprinted on the swarm orientation and can cause the formation of filaments and patches (see panels (d-g) in Fig 6 and S5 Movie).
We quantify the success in avoiding turbulent regions by means of the indication number Eq (1) and define the fraction: When the agents all avoid turbulence, PðjCjÞ ¼ 0, when agents are initially randomly distributed and the active agents all reside inside areas of high turbulence, P ¼ 1=N ðjCj passive Þ.
The characteristic regions where jCj > jCj passive;N are indicated by the black boundary in Fig 6. We see that these regions for the cavity flow evolve essentially into two substructures, an elliptic central region and a long and narrow hose-like structure generated by the clockwise main current.
Mainly influenced by these flow structures, there appear essentially two classes of transients. First, when counter-current swimming is switched off or if the orientation against the flow is weak (ε U = 0.005 rad s −1 ), the number of agents in turbulent regions decreases with increased ε C and the transition out of turbulence takes place essentially in a short initial time interval (see Fig 6 panels (a)). Second, when counter-current swimming is sufficiently strong, the jellyfish effectively stand still in the water or achieve a net velocity against the flow. In that situation, they stay in a surrounding turbulence structure for a prolonged period of time and can group along this structure into patches or long filaments (see Fig 6 panels (d-g), S5 Movie). Nevertheless, the swarm can eventually encounter turbulent regions in which it stays for some time before it starts avoiding these regions (see Fig 6 panel (f)). When this happens, we observe intermittent spikes in PðjCjÞ that decrease in amplitude before the swarm has adjusted to the turbulence structures (see Fig 6 panel (b)).
Since spikes are isolated and transient events in time, we first omit the initial 50 s and calculate the maximal and minimal fractionsPðjCjÞ andPðjCjÞ. Then, we take into account only values PðjCjÞ < 0:9PðjCjÞ þ 0:1PðjCjÞ to obtain the time average P min ðjCjÞ (see Fig 6 panel  (c)). The resulting curves show that the baseline avoidance of turbulence is essentially independent of ε U,0 .
Additionally, we characterize the pattern formation by ΔHex 12 (X) = hhHex 12 (X) − Hex passive,12 (X)ii t and the difference of standard deviation ΔS = hhS(X) − S passive (X)ii t where Here, σ x and σ y are the standard deviations of the agents positions in x and y-direction respectively at time t. At t = 0, S(X) = 1 and stays approximately constant for the passive swarm

PLOS ONE
while it increases or decreases according to the evolution of the active swarm. We simulate the dynamics for 900 s and take into account values in the last 300 s for time averaging. In contrast to the analysis of regular structural order, we have chosen 12 nearest neighbors for averaging of the phase factor exp(i6δ) in Eq (18), because it allows us to better resolve the effect of filamentation.
We observe that the structure formation process heavily depends on the combination of ε C;0 and ε U,0 . When swimming against the flow is switched off or if it is weak, the jellyfish follow scattered trajectories and mostly get trapped inside the corners of the domain, causing a large spread and a relatively high order of the swarm due to compression at the walls (Fig 7  black and blue curves and Fig 8 panels (i, j, m, n), S6-S9 Movies).
With increased ε C;0 the swarm transitions into a single patch in the center (Fig 8 panels (h,  k, l, o, p), S10-S14 Movies) which has a minimal spread and a minimal order. The latter is a consequence of the constant regrouping due to agents that bump into the cluster and advection. On the contrary, when the agents orient more against the flow, they are able to group into long lasting and coexisting coherent structures, namely filaments (Fig 8a-8g, S15-S21 Movies) and rings (Fig 8a, 8b, 8e and 8f, S15-S20 Movies). Most prominent is the transition for the optimal, ε U,0 = 0.16 rad s −1 . For this setting, the swarm undergoes a transition from ring patterns over filaments towards patches. This transition is indeed signified by two jumps in the spread ΔS and an increase of local order before the system settles into the disordered patch (see green lines in Fig 7). In contrast, for ε U,0 = 1 rad s −1 , the swarm transitions later from ring patterns into filaments. Based on this investigation, we chose ε C ¼ 0:08 rad s −1 (ε U ¼ 0:16 rad s −1 ) as an optimal pair of coupling values (red diamond in Fig 7, see also Fig 9). It gives rise to a dynamics that sits "in between" scattered states and patches and thus, ensures a relatively rich behavior of the simulation model as already small variations of ε C will lead to a significant change in the swarming patterns. We again favor such a behavior as it mimics the frequently encountered phenomenon of criticality in biological systems [138,139].

Foraging
Our final goal is to capture the principal behavior of foraging based on the full dynamics Eqs (3)-(11) and (13)- (16). We simulate the swarm in a double-gyre flow in which we initialize the agents in the right gyre at x > 0 and the prey in the left gyre at x < 0 (see Fig 10 panel (a)). Thus, in this setting the jellyfish are dynamically separated from prey by a strong central current at x = 0 [35]. We start our model with the previously obtained optimal parameter setting

PLOS ONE
of Table 2 for counter-current swimming and turbulence avoidance. Overall, each model instance runs for 900 s.
Response mechanisms. The turbulence avoidance causes the swarm to form a ring structure in the right gyre of the flow (see Fig 10 panel (b) and Fig 11 panel (a)). We check first the impact of subsequent combinations of response mechanisms on the ability to search for prey: A Higher pulsation frequency oðAÞ and a larger angular diffusion DðAÞ.
B Reduction of turbulence avoidance by means of a response ε C ðAÞ. We set the slope parameter for ε C to f C ¼ 0:02, allowing for an almost immediate reduction when activity increases.
C Bi-linear dependence of propulsion velocity VðU j þ A j Þ.
D The full response dynamics at optimal parameters of Table 2. (This search dynamics is depicted in Fig 10 and can be observed in S22 Movie).
Prey search by jellyfish is variable but there exists experimental evidence [35, 105] and observational experience [65] that can be condensed to the following two performance criteria for our tank simulation: When jellyfish is searching successfully for prey,

PLOS ONE
• The fraction of jellyfish in the right (oligotropic) gyre is minimal, meaning that the majority of individuals manages to cross into the left domain at x < 0 or that individuals group in regions where prey is found.

PLOS ONE
• The jellyfish spread evenly in the left gyre once they have crossed the barrier. This mimics the free floating state when individuals shuffle water to the oral arms.
We quantify the crossing performance by means of the fraction where we indicate by index 0 that only jellyfish in the right gyre (x > 0) contribute (x ¼ 0). Accordingly, a fraction PðxÞ ¼ 1 indicates that all jellyfish and all passive agents are found in the right gyre and it is zero, when all jellyfish have crossed the flow barrier. To quantify the spread of the swarm after the crossing of the barrier, we calculate S 0 (X) according to Eq (22). The index 0 indicates that for this averaging only jellyfish at x < 0 are taken into account. In this analysis we do not take into account the length of the transition time needed for the swarm to switch basins. However, we found that the transiton, if it takes place, is completed in the end of the simulation time. While frequency and diffusion responses (A) do not result in significant searching (Fig 11  panel (a), S23 Movie), already enabling weakening of the turbulence avoidance (B) allows some agents to cross the flow barrier and to spread partly in the left gyre (Fig 11 panel (b), S24 Movie). A velocity response to activity increase (C) leads to a further increase in the number of crossings. However, jellyfish get trapped in the corners of the tank because direction is no

PLOS ONE
longer chosen based upon turbulence avoidance and higher velocities increase the chance of bumping into a wall where the movement is constrained. (Fig 11 panel (c,d,e)), S25 and S26 Movies).
Interplay of directed search and counter-current swimming. Next, we also allow for an orientation into the direction of the prey. We find that this additional orientation (δ F ) leads to complex decision making in our jellyfish model based on an interplay of directed prey search and counter-current swimming. We start our analysis by limiting the set of prey angular-coupling constant and response parameters ½ε F;0 ; f C ; f F ; f U ; f g � to symmetric switching from turbulence avoidance to prey searching (f F ¼ f C ). And, we fix the ignorance parameter f g = 100, which means that even after prolonged periods of preying, resulting in high activity A, the jellyfish remain "greedy" for food. This leaves us with ε F,0 , f F and f U .
The ability of the swarm to cross the flow barrier is found to depend mainly on the angular coupling constant ε F,0 and velocity decoupling parameter f U . On the contrary, crossing performance is largely independent of f F (see Fig 12). In particular, we see that for f U = 0.05F 0 (86% weakening of counter current coupling at F j = 0.3 F 0 ) more jellyfish travel to the left gyre ( Fig  12 panel (a)) than for f U = 0.01F 0 (97% weakening of counter current coupling at F j = 0.3 F 0 ). Interestingly, this difference occurs due to the emergence of secondary clusters in the right part of the tank (see for instance panels (b,c,f) in Fig 13). These clusters emerge due to a longtime equilibrium of orientations against the flow (δ U ) and towards prey (δ F ). The upper cluster is less stable than the lower one because jellyfish in the upper right part of the flow are the last ones to become activated. As a consequence, they remain too slow for the flow in that region and experience a slow downstream transport. Ultimately, they come close to the upper central region of the tank where they feel the prey and cross the barrier.
On the contrary, agents of the lower cluster experience a more diffuse flow and thus, are able to counteract advection. Ultimately they are able to approach the flow barrier where an increase in activity sparks a shift in the coupling parameters such that they cross over as well. However, those jellyfish also experience a small leakage of prey from left to right domain, which is transported by the lower boundary current in which they swim. This, on the one hand causes stronger bell strokes. On the other hand, the jellyfish start to orient partly towards the ambient prey. But, since the prey gradient and the direction of the nutritious flow are mostly not aligned, many jellyfish get dragged away. Some of them find a new equilibrium position further away from the flow barrier while others get transported with the current to the upper part of the barrier where they cross. Overall, some of the jellyfish remain inside the nutritious current, mainly on the lower right side of the tank.

PLOS ONE
In fact, coupling values ε F,0 > 0.16 rad s −1 make the escape from the lower secondary cluster more difficult while at the same time causing the swarm to have a small spread in the left domain (Fig 12 panel (b)). On the contrary, for ε F,0 < 0.16 rad s −1 , the jellyfish do not orient sufficiently towards prey and get trapped in the corners of the tank. This leaves us with a narrow parameter region of optimal ε F,0 (grey regions in Fig 12) out of which we choose ε F,0 = 0.16 rad s −1 .
Finally, we investigate how the response to prey depends on the velocity decoupling parameter f U . For this we consider just ε F,0 = 0.16 rad s −1 . Here, we find that because the value of f U determines the instantaneous coupling strength ε U (F), it simultaneously influences the fraction hPi of remaining agents in the right gyre and the swarm spreading in the left gyre (Fig 14  panels (a-d)). We see that for largely persistent counter-current swimming almost all agents cross the flow barrier (f U = 0.5F 0 , 38% weakening of counter current coupling at prey concentration F j = 0.3F 0 , Fig 14 panel (d), S27 Movie). However they do not spread in the left domain but form a ring pattern. On the contrary, for earlier loss of counter-current orientation (f U = [0.025, 0.075, 0.1] F 0 , [92,80,75]% weakening of counter current coupling at F j = 0.3F 0 , Fig 14  panel (a,b,c), S28-S30 Movies), secondary clusters remain in the right domain (Fig 14 panels  (a-c)) but the ring patter in the left gyre is less pronounced. Overall, this leads us to the conclusion that we consider f U = 0.05F 0 as an optimal parameter as it allows almost complete crossing of the swarm and reasonable spreading of the swarm in the left domain of the tank (see Fig  10).

Conclusion
In this paper we have proposed a paradigmatic model for jellyfish swarming based on active Brownian particles. The model can be readily incorporated in large-scale simulation models of the ocean and simulations of tank experiments. It provides a measurement paradigm to understand mechanisms of large scale structure formation based on agent responses, agent physiology and local interaction mechanisms. In its present form, Eqs (3)-(11) and (13)-(16), the model is already able to reproduce paradigmatic jellyfish behavior reported in the literature. In particular, we have been able to simulate the counter-current swimming of R. nomadica [23]. Using the kinematic tracer of absolute vorticity, we suggested a mechanism for the formation of long-lasting swarm patterns. Finally, we simulated the qualitatively reported behavior of jellyfish in a tank when prey is present [35, 75,105].
Here, we will discuss several aspects of swarming that have to be investigated in a future model. We largely dismissed jellyfish physiology of swimming and agent-agent interactions. The physiology of a jellyfish comprises three aspects. First, we expect that major parts of communication among jellyfish can be attributed to the turbulent vortex rings which are send into the water with a certain angular dependence Q V (x k , x j , θ k , θ j ) (red arrow in Fig 1 panel (a)).
Second, jellyfish possesses specific sensing capabilities. It has been reported that jellyfish occasionally bump into each other "head on" while at the same time being able to sense already faint pressure signals at a distance [65]. This suggests that jellyfish perceive their surrounding depending on the distribution of rophalia on the surface of the bell, introducing an angular  Table 2. Color indicates the activity A. Black dots indicate the swarm in the end of the simulation. Panel (e) depicts the time averages of the agents fraction hhPðxÞii t at x > 0 (black squares, left ordinate) and of the spreading hhS 0 (X)ii t when x < 0 (blue dots, right ordinate). https://doi.org/10.1371/journal.pone.0288378.g014

PLOS ONE
dependence of sensitivity to perturbations. In the given framework, we suggest to introduce a function Q S (x k , x j , θ k , θ j ) that corresponds to this distribution. Accordingly, the resulting physiological effect on agent-agent coupling, can be incorporated by a product function Q(x k , x j , θ k , θ j ) = Q S (x k , x j , θ k , θ j )� Q V (x k , x j , θ k , θ j ). The function Q(.) then mimics the previously envisaged angular dependence of coupling in attraction and repulsion Eq (16) and in the functions H j (.) and G j ð:Þ, Eq (2).
Third, one needs to specify how the dynamics of the internal neuronal network reacts to external stimuli. This effect can be taken into account by the phase coupling function H j (.) of a medusa. It describes the response of the network phase, φ, Eq (2) to external perturbations and has found widespread application in research on biological oscillators [46,70,103,115,[140][141][142]. As a consequence of the directional dependence of agent-agent coupling, given by Q(.), the phase response and the sensitivity response have to be disentangled simultaneously. Here, tracking techniques and phase-dynamics reconstruction have to be employed [143][144][145].
Currently, we have considered only a simplified fluid-like soft-core repulsion that avoids any attractive effects. Here, further investigation of the radial dependence of the flow field around a medusa are needed. Beyond the angular dependence of communication, given by Q (.), jellyfish might also be subject to inevitable effects of synchronization among each other if it fits their energy household [124,146]. This effect would involve the oscillatory phase of the bell, φ and further modifies the attraction and repulsion in Eq (16).
Generally, synchronization effects can be subtle because biological oscillators tend to only partially synchronize. From that perspective it is an intriguing open question to which extend synchronization of bell pulsation is present in jellyfish aggregations and whether it leads to changes in swimming performance of jellyfish in blooms. For example, in our model the almost complete absence of a mechanism that could lead to synchrony of the bell oscillations causes significant performance losses of jellyfish clusters due to jamming (for instance S11 Movie). In fact, the only mechanisms through which the bell oscillations can couple in the model is present when agents search for prey. We consider this effect to be insignificant as it is caused by the interplay of prey concentration, repulsive interaction, activity dynamics and frequency response. On the contrary, bells stay effectively asynchronous such that densely packed agents push each other out of the way until they have sorted into groups where pushing becomes rare. This mechanism is know as social sorting [63]. Indeed, this effect might be partly visible in S10, S12 and S14 Movies where agents form highly dense clusters. In our given model, walls and agents have the same effect on the positional dynamics, depending on the interaction distance R i . We have seen during our numerical analysis that the choice of R i significantly influences the foraging of jellyfish in the secondary clusters at the walls. Thus, a further open question is the nature of the wall interactions for jellyfish.
We have completely ignored perturbations in the fluid due to the presence of agents. Our motivation here is that our model will serve as a parameterization of a sub-grid-scale process [147][148][149]. As such, the local fluid perturbations of single jellyfish become irrelevant. Therefore, we take the local perturbations into account only by the angular noise L j . One of the two essential parameters of this noise, D 0 , was chosen and thus, introduced a time scale to which the coupling constants of orientation relate. Moreover, we use a normalized prey concentration by introducing the maximum level of prey F 0 . To make the simulation more accurate, the actual statistical properties of the angular noise need to be investigated from passive observations and the prey concentration needs to be normalized according to an empiric F 0 . Furthermore, in a large-scale ocean simulation, the fluid-agent coupling will require a high-order interpolation scheme in time and space to provide accurate Lagrangian field values. In turn, the agent-fluid interaction shall rely on averaging over agents in a single cell.
The decision making of jellyfish in our model is realized by response functions Rð:Þ that differ only by two parameters and are otherwise similar in shape (see panel (c) Fig 3). We have employed this concept of response functions because of its great success in the description of cell membrane responses in computational neuroscience [106,107]. In contrast to those systems where the membrane response was measurable directly, investigation of the response parameters will have to rely on indirect observations and data fitting methods such as machine learning [72].
We have simulated a two-dimensional swarming dynamics. Thus, in its present form, the model Eqs (3)-(11) and (13)-(16) can be used to capture the vertically averaged dynamics of a swarm and the dynamics of jellyfish in periods of reduced vertical migration. However, in reality, jellyfish swimming is influenced by the day night-cycle [88,150], its prey and environmental inputs such as temperature and salinity [13] which cause vertical migration. Augmentation of the model dynamics to three dimensions is relatively straight forward and requires a second directional angle.
We have chosen an over-damped positional dynamics Eq (14) to model the swimming. By this we have in mind, that on larger scale, in the sea, the short time and length scales associated with inertia are not relevant. However, this assumption gives rise to two effects in our work. First, we observe regaining of mobility once the density of individuals becomes large enough (see Fig 4). Second, the phase in the velocity oscillation function β(φ) Eq (15) has only minimal effect on the positional dynamics as the time-average of velocity is completely determined by the shape parameter J. A future model might consider also inertia for jellyfish, including their velocity dynamics [42,53]. However, such a model requires twice as much initial conditions (for the positions x j (0) and velocities v j (0)), imposing even greater demands on any datadriven set of initial conditions.
The turbulence avoidance and the prey searching represent a type of first-passage problem in which the supremum of the first passage times of single agents could be used to define a first passage time for the swarm to enter the target region [151]. However, we would like to point out that agent trajectories are not independent of each other and that the boundary of the target region can be expected to have a fractal, time-dependent structure according to the Reynolds number.
Finally, we have investigated just a finite time horizon of maximal 15 minutes due to numerical constraints. We think that in future analysis longer observation times are needed to capture long-term effects of swarming.